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The classical limit of the Wigner- Weyl representation is used to approximate products of bound- 
• continuum matrix elements that are fundamental to many coherent control computations. The range 

of utility of the method is quantified through an examination of model problems, single-channel Na2 
i— i, dissociation and multi-arrangement channel photodissociation of CtblBr. Very good agreement 

with the exact quantum results is found for a wide range of system parameters. 

I. INTRODUCTION 

<D ■ 

Recent experimental and theoretical developments show that coherent control, i.e. the control of atomic and 
, molecular dynamics via quantum interference, can be successfully applied to the wide variety of systems 0. However, 
O ' available quantum methods limit theoretical studies to scattering and photodissociation of small molecules. The 
extension of theoretical studies to larger complex chemical systems must rely on new developments in quantum or 
£>v semiclassical techniques. Primary amongst these is the application of the Initial Value Representation semiclassical 
' approach that is described elsewhere However, the ultimate utility of such techniques requires the development 

i ^~f , of numerical tools to speed convergence of the initial value integrals with highly oscillatory integrands. 

In this paper we consider an approach to computing interference contributions in coherent control in which such 
oscillations do not occur. Specifically, we focus attention on using the classical limit of Wigner Phase Space Methods0, 
0, where the desired transition matrix elements products are written in the Wigner- Weyl representation and the 
quantum expressions are replaced by their classical counterparts. This approach has been applied in the past to a 
number of simplerproblemsj^, Q where transitions from one initial bound state were considered, e.g., from a single 
ly-^ vibrational state 6] or from a Gaussian state 0, @ • In these cases only absolute values squared of matrix elements 
were sought. By contrast, in this paper we use a result [lfj on the classical limit of the nonstationary Liouville 
(*C) ' eigenstates to obtain the classical limit of the desired product of transition dipole matrix elements, a quantity which, 
in general, includes phase information. For example, we demonstrate that the method gives good results for complex 
valued matrix elements that arise in the multi-arrangement photodissociation of CH2iBr to CH2I + Br or CH2Br 
+ I. The method has also been previously applied in the guise of the linearized approximation to the Initial Value 
* 55 ' Representation with varying degrees of success. 

This paper is organized as follows. The classical limit of the Wigner phase space method is described in Section 
II. Applications of the method to simple systems involving one arrangement channel are described in Section III. 
Specifically, we consider Franck-Condon transitions for the model case of excitations from a harmonic oscillator 
potential to a linear potential and to transitions on realistic Na2 potential energy surfaces. Section IV discusses 
applications to collinear photodissociation of CH2lBr, where the transitions probabilities are complex because the 
relevant operators are no longer Hermitian. Section V provides a summary of results. 
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II. METHOD 



In various coherent control scenarios (e.g., bichromatic control 12] or coherent control via pulse sequencing [l3|1 
control is dictated (here written for the case of non-rotating diatomics) by interference terms of the form 

k k 
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Here \xj) = (J-fi\j), j — n,m and fj.fi is the electronic transition dipole moment between the upper and lower electronic 
states. The states | j) of energy Ej, are bound vibrational states on the lower potential energy surface, and | E, k, r~ ) 
are continuum nuclear eigenstates of energy E on an excited electronic surface, corresponding to product in channel 
r with quantum numbers k. Terms like those in Eq. arise when one calculates the probability P(E) of transition 
from a coherently prepared bound superposition state 

w = 5>i*>, ( 2 ) 

i 

to the final continuum states at energy E in arrangement r |T3 |. That is, these terms correspond to the interference 
terms in computing the probability P(E) oc J^k I ^ r_ Im/jIV') | 2 - Since coherent control relies upon quantum 
interference to obtain control over molecular outcomes, the evaluation of terms like those in Eq. Q are vital to 
computational control studies. 

In the simplest case, when n = m and only one product arrangement channel is open, Eq. is proportional to 
the photodissociation probability for the transition from an initial bound state |n) at energy E n to the continuum of 
final states at energy E: 

*£ n (E)=J2\( E ^r-\» fi \n}\ 2 . (3) 
k 

To utilize Eq. we rewrite it as: 

a^} m {E) = Tr[(]T | E, k, r~ ){ E, k, r" \(\ Xm ) (Xn\)}, (4) 
k 

where "Tr" denotes the trace. The term J2\t I E, k, r~ ) ( -E, K, r _ | can then be written as 

^|£,k,r)(£,k,r" | = R r S{E - H) (5) 
k 

where i? r projects onto product arrangement r and H is the Hamiltonian of the excited electronic state [T^J. When 
summed over r (or when there is only one product arrangement channel) this equation reduces to the familiar 
expression: 

J2\E,k,r-)(E,k,r-\=S(E-H). (6) 

k,r 

Using Eq. J5J, Eq. (0} can then be rewritten as 



v£UE) - Tr[R r S(E - H)}(\ Xm )(Xr, 



(7) 



The Wigner transform Oiy of any operator O is defined as 

r+oo 



1 f + 

Ow(p,q) = -r / dv e 2i P- v / R (q- v|0|q + v), 



(8) 



where p is the momentum conjugate to coordinate q. Taking the Wigner transform of Eq. J7|) and using Tr(AB) = 
Ti(A w B w ) gives 



(E) = Tv [R r S(E~H)] w p 



X 

n.m 



(9) 



where p^ m (p, q) = [|Xn)(Xm|]w is the Wigner transform of \Xn)(Xm\- Neglecting the coordinate dependence of fj.fi 
in accord with the Franck-Condon approximation gives 



(10) 



where p n ,m = Pn,m(p, q) = [| n )( m Equation I|1U|I is exact. In Section III of this paper we consider applications 
of approximations to Eq. I|1U[) to cases with a single product arrangement channel. Section IV considers the case of 
multiple chemical products. 
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III. SINGLE PRODUCT ARRANGEMENT CHANNEL 

In the case where there is only a single product arrangement channel, Eq. 1|10[) becomes 

cr n , m ( E ) = p}i Tt [ S ( e - H)wPn,m] = M/j J dp dq 5(E - H) w p n , m (p, q) . (11) 

Pn,m(p,q) = [|n)(m|] w = ^ J dv exp(2ip- v/7i)(q- v|n)(m|q + v) , (12) 

so that a n ,m{E) appears as the overlap of two phase space densities, one corresponding to the density on the lower 
surface [p n . m (p, q)] and one to the continuum density 6(E — -ff)w on the upper excited surface. Neither of these 
densities need be classical since they can be negative or complex. 

The classical Wigner approximation to a n , m (E) is obtained by taking the classical limit (h — > 0) of Eq. i|ll[l . 
where the classical limit of 5{E — -ff)w arises by either expanding the density of states in powers of Ti, or by using 
the statistical operator dH P = exp(-f3H), f3 = 1/kT, or via an exponentiated h expansion [T^. or by expanding 
8(E — H)w around the identity operator / times the classical Hamiltonian, H(p, q) ■ / ||. The lowest order term 
in this expansion is S(E — H(p,q)), where H(p,q) is the classical Hamiltonian associated with the upper potential 
energy surface V(q). For a particle of reduced mass p in one dimension, which we focus on in this section: 

H(p,q) = £+V(q) ) (13) 

so that the lowest order classical Wigner phase space method approximation to <T n ,m(E) [denoted o c n rn (E)] is 

< m (£) = fx% [ dpdq 6(E - H(p,q))p c ntm (p,q) (14) 



f ,1 A ~ l(P' E )) c I N 2 
j / dpdq — Pn,mXPi1) = M/i 



(*) 



q=q(p,E) 



dp 



Pn (q(p,E),p) 



(f) 



q=q(p,E) 



(15) 



where E — p 2 /2/x — V(q(p, E)) = 0, and where Pn, m (P> l) ^ s * ne classical limit of p n ,m ■ 

The form of p c n m depends on whether the system is integrable or non-integrable j 1 0| . For integrable systems [Eq. 
JHJ, Eq. CHI), Ref. [13] 



1 

2^ 



S(I(p, q) - J n m ) exp[i(n - m)9(p, q)), 



(16) 



where [I(p, q), 9{p, q)] are action-angle variables, I n m = (I n + I m )/2, and I n is the semiclassical action associated with 
state |n) (i.e., I n = (n + "/)h, where 7 is the Maslov index). A similar analytic expression is not possible for chaotic 
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systems 

In the cases studied below we focus on transitions from low lying vibrational states of diatomics. In this case we 
can approximate the potential by an harmonic oscillator of frequency to and approximate p c n m (p, q) by the harmonic 
case. However, in the harmonic case the classical p c nrn {p,q) can be chosen as equal to p n ,m{p,q) for the quantum 
harmonic oscillator 16]. Substituting expressions for the harmonic oscillator states 



(q\n) = (^j (2»n!V^) 
into the Wigner transform of Eq. I|12|) , we obtain 

(-1)" 



1/2 HjJlf q] e- 



(17) 



P C n ,m(Pi<l) = Pn,m(P,q) = 



ttTi 



2™m! 
n\2 m 



1/2 



fipu 



[ip — pujq\ exp I I L, 



4I(p,q) 



(18) 
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where n > m, is the generalized Laguerre polynomial, and the classical action I(j>, q) for the harmonic oscillator is 

H(p,q) 1 ( p 2 2 \ 
I(P, 1) = = o + ■ (19) 

UJ 2 \pU! J 

For the case of n — to, Eq. 118(1 takes the well-known form 

(-1)" / 2H(p,g) \ ( AH(p,q) \ 
Pn,mip, q) = —t- exp — — L n — . (20) 



IV. RESULTS: SINGLE ARRANGEMENT CHANNEL 



A. Model Potentials 



We first test the utility of this approximation on a simple standard model: excitation from an harmonic oscillator 
initial state potential to a linear, repulsive excited state potential V(q) — —j3q + E , with E arbitrary and > 0. 
This model was previously examined for n = to = in Ref. Q, and can be used to approximate transitions to an 
arbitrary potential if is taken as the slope of the upper potential energy surface at the peak of the initial state 
(n = m = 0) wavefunction. For simplicity, we set the pfi to unity. 8{E — H)w is known |6( for the linear excited 
state potential, so that the exact (J n ,miE) is given by 

a n , m (E) = 2ir dpdq5(E - H) w p n ^ m (p,q) 



-(2p 1 ^){q~q(p,E))\ Pn ,mip,q) 



(21) 



where Ai is the Airy function. 

By contrast, the classical result [Eq. 1(15(1 ] for this model is 



< m (£) = J dpdqsfE-^+Pq-EhjpZ 



= i J dpdqS(q - q(p, E)) p^ m (p, l) = \J dp p c n , m ip, q(P, E)), (22) 

where q(p, E) = (E - E)/f3 +p 2 /2p(3. 

In Fig. 1 we compare the exact quantum and classical Wigner results for the highly quantum case of p — m e —1 
a.u., where m e is the mass of electron. Note that <7 n>m , n ^ m is real for the one-dimensional case with real V{q) 
since the integral over the imaginary part is odd in the momentum variable. In Figs. l(a)-(b), a n ^ m for n = m = 
and n = m = 4 are shown as functions of energy E — Eo for parameters given in Ref. 0; (3 = 6 a.u., uj — 2 a.u. Our 
results for n — agree very well with those in Ref. |fg, and with the n = m = 1 results for = 6 a.u., u> = 1 a.u. 
(not shown). However, as is evident from Fig. 1, the accuracy of the classical Wigner approximation deteriorates 
extremely rapidly with increasing n; results are very poor even for n = 4. The same behavior is observed for cases 
where n 7^ to, shown in Fig. l(c)-(d) for (n = 1, m = 0) and (n — 5, m — 4). This is because for small n, p n , n 
is smooth and broad, and the transition integral averages over many of Airy function oscillations. By contrast, for 
large n, p n ^ n is highly oscillatory and the initial state probes fine details of the Airy functions which are absent in the 
classical approximation. As is evident from Fig. 1(b), the classical Wigner results for n — to are negative at some 
points. This is impossible physically and indicative of errors in the approximation. In the n ^ to case, where negative 
values are possible [Fig. l(c)-(d)], the positions of maxima in the semiclassical approximations are shifted to lower 
energies, a feature previously explained for the n = m case 

Results presented in Fig. 1 were obtained for p = 1 a.u. which is highly non-classical system. The dependence of 
the classical-quantum agreement on the oscillator frequency lo and on the reduced mass p are explored, for two cases 
with to = n, in Fig. 2. Agreement is seen to improve dramatically as to decreases and somewhat as p decreases. The 
latter result is surprising, motivating the analysis described below. 
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Characteristic s n factor. The validity of the classical Wigner approximation for cases where n = m and for the 
excited state linear potential has been examined in Ref. There, the utility of the Wigner phase space approximation 
for an initial Gaussian wavefunction was shown to depend upon two parameters, s and A: 



1/3 



where A sets the scale for the width of the excited state wavefunction's oscillations near the turning point (other 
oscillations have shorter wavelength), and A is the width of the initial Gaussian ground state: 

^^-^m^ (q - qo)2,2A2 - ( 24 ) 

Specifically, Eckhardt and coworkers have shown that the smaller the s, the more accurate the classical approxima- 
tion, consistent with the fact that smaller s means more Airy oscillations over the width A. However, the parameter s, 
as defined in Eq. i|23|) . is not as useful in our case since we consider assorted ground state vibrational wave functions, 
and not just simple Gaussians. Specifically, s was unable to predict the correct dependence of the classical-quantum 
agreement on \i or to. This is because for any oscillator, the width of the vibrational state A depends on both \x and 
u>, whereas Eq. (|23|) only allows for a dependence on \i via A. 

To generalize the s expression to higher vibrational states, for cases with m = n, we compare the known expression 
for the v — vibrational level 

i'= [-) e 2 , a = -p, (25) 
to Eq. I]24p. and obtain A = -7= = y^j- The s parameter can therefore be written as 

S= A = \-W ■ (6) 
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Further, to account for the oscillatory character of the nth wavefunction on the lower potential surface, A is replaced 
by A„, where A„ equals the width at half- maximum of a single oscillation of the ground state wavefunction. This 
is given by A„ = l/2(n + 1), where I = 2^j2{n + T/2jKJJuv is the overall width of the ground state wavefunction, 
estimated as the distance between the two classical turning points. This gives the new parameter s n : 



A / (n + 1) 2 fn , . fntT . 

- • or s„ = W- s for large n. (27) 



" A„ y2(n + l/2) ' " V2 

It follows from Eq. (|27|) that the smaller the s is for n = m = 0, the better the classical approximation for higher-lying 
levels. However, s n will always increase with increasing n due to the decreasing wavelength of the bound wavefunction. 
Further, decreasing the vibrational frequency of the lower electronic state u> leads to wider vibrational states (and 
larger A in Eq. (|23[1 for the Gaussian ground state) and, therefore, to the smaller s and better agreement with 
quantum result. This explains the results shown in Fig. 2(a)- (b), where the corresponding values of s for two u> values 
are 0.62 (s4 = 1.03) and 0.20 (54 = 0.3), respectively. Alternatively, note that for smaller ui, the density p n , m becomes 
wider in q and the overlap integral in Eq. Ijlll) averages over more Airy function oscillations. Increasing the slope of 
the upper potential also leads to better agreement between the classical Wigner and quantum results in accord with 
Eq. 1261) . reflecting the fact that the Airy function oscillates more with increasing f3 [Eq. I|21l) ]. Further, Eq. 1|26|) 
and Eq. (|27l) quantify the dependence on reduced mass [i that is evident in Fig. 2(c)-(d) where corresponding values 
of parameter s for the two increasing values of u are 0.62 and 1.33, respectively. This is because, while A in Eq. ll'-il) 
decreases with increasing /i, the width of the initial state A does so as well. 

Characteristic s n ,m factor. In the case of n 7^ m, developing the corresponding parameter s n ^ m is more difficult. 
However, in accord with s„, the quantity s n ^ m is expected to be the ratio of a characteristic width on the excited 
state divided by a width on the ground state. To obtain s„ iTO we make the following observations: (1) The width 
l m = 2y/2(m + l/2)h/fuv of the product (f> n 4>m can be defined by the width of the narrower of the two states m, with 
m < n, where <f) n is given by Eq. (jT7|) . since (f> m — > outside that interval; (2) the total number of relevant nodes of 
the product 4> n 4> m is n' + m' , where n' and m! are the number of zeroes of the nth and mth states within the interval 
l m . Further, the quantity m! — m since l m is the overall width of the mth harmonic oscillator state and n' can be 
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estimated as l m /A n , where A„ is the characteristic width of cj> n ; 



a jn _ 2^2(n+l/2)ft//^ 

™ Ti+l rc + 1 1 j 

Hence, the number TV of oscillations of 4> n 4>m is given by 

N=^+m + l. (29) 
Thus, the characteristic width at half-maximum A n m of the oscillations of the product </> n </> m is 



_ l m /2 _ l m /2 _ y/2(m+l/2)h/pv 

L±n.m — 77 — ~j 7 — ; 

^-+m+l (n + 1) f*jiE + m + 1 



The parameter s„. m is then 



A y^(n + l)+m+l 
S "' m -A n , m " (m + 1) Sm - (31) 



Substituting 

A m + 1 s 



(32) 



we obtain the pleasing result that 

Sn,m — 1/2 \S n + 5 m ) . (33) 



As expected, for n = m, s rhn — s n . 

The behavior of the results in Fig. 1 can now be quantified in terms of s n and s n>m . In particular, we obtain 
the following values of s n and s n ,m for the results presented in the figures: So = s = 0.62, S4 = 1-03, si : o = 0.67, 
and S5 ; 4 = 1.08. One can see that s is already larger than 1 (fast diverging series; see Table 1 in J8j) for n = 4. 
Additional computations show that s^o = 0.83 and the approximation works reasonably well in this case. Similar 
good agreement has been obtained for three different pairs of n and m which are characterized by s n>m w 0.83, e.g., 
(n,m) — (4,0), (3,1), and (2,2). However, in all cases where s„ iTO > 1 (e.g., sal. 12 for (n,m) = (8,2), (7,3), and (5,5)), 
the agreement is poor (not shown). 



B. Molecules 



To apply this to realistic systems, consider first a linear potential model of the dissociation of the H2 molecule. 
In this case \x is increased and the ground state uj is decreased relative to the system studied above. In particular, 
/j = 918.7 a.u. and u> = uj e — 4395.2 cm" 1 17]. o results for four pairs of n and m are shown in Fig. 3 for (3 = 
6: n = m = 4, n = m = 20, n = 5,m = 4, and n = 21, m = 12. The results clearly show much better agreement 
between exact Franck-Condon and classical results than does the case studied above. Indeed, the difference between 
the classical and exact quantum results becomes visible only at high-lying levels, n = m = 20. Furthermore, this 
difference becomes practically indiscernible for the even heavier molecule Na 2 , fj. — 20953.9 a.u. uj — uj^a 2 = 158.91 
cm -1 0], where virtually perfect agreement between exact and approximate results is obtained even for n = 20 (not 
shown). Using Eq. (|26J) we obtain the corresponding values of parameter s: S[5] = 0.62 (i.e., the parameters are 
those of Ref. [fj and used in Section IV. A above), sh 2 = 0.19, and SNa 2 = 0.06, in accord with the computed results. 
Similarly, the excellent agreement obtained for the high n and m values for model H2 and Na2 is consistent with the 
values of s n and s n ,m- For example, S20= 0.197 for Na2 and 521,12 = 0.179. In the case of the H2 model S20 = 0.62 
and we expect to sec the difference between the classical and exact quantum results on a par with the one we have 
seen for the case in Section IV. A where sq = s was also 0.62. This is indeed the case, as seen from a comparison of 
Figs. 1 and 3. Note also that S2i,i2 = 0.57 in the case of the H2 model which leads to worse agreement than in the 
Na2 case (s^i 12 = 0.18) but much better than in the case in Section IV. A (S2L12 = 1-85) (not shown here). 

To see the origin of this behavior, consider cuts through the ground and excited phase space distributions at fixed 
values p = and E = Eq, as shown in Fig. 4 for three different systems: the model of Ref. |6J, H2 and Na2. While 
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FIG. 3: Comparison between the classical Wigner (solid line) and exact quantum (dashed line) results for a n , m (E — Eq) (n = m) 
with /i = /xh 2 and o> = <^h 2 — 4395.2 cm -1 : (a) n — m — 4; (b) n — m — 20; (c) n = 5, m = 4; (d) n = 21, 12. 



the quantum results correspond to taking the overlap between these two phase space densities, the classical limit 
corresponds to using the value of the initial state phase space distribution at the coordinate q(p, E) (the vertical line 
in Fig. 4 corresponds to the value q(p = 0, E = Eq) = for p — 0, E = Eo)- Clearly, this approximation improves 
with decreasing s. 

In an attempt to further improve these results, we considered the higher order quantum corrections (see Eq. (4.8), 
Ref. 0) to the classical Wigner result, which consist of additional terms in an expansion of the density of states 
S(E — H)w m powers of %. Our calculations using this expansion showed that although these corrections sometimes 
work well for the low- lying vibrational levels n < 2, they led to much poorer agreement with the exact results for 
higher n, as the higher order terms became dominant, a result also noted in Ref. 8]. As is evident from Fig. 5, 
where the quantum corrected results are shown as dot-slashed curves, the quantum corrections are practically of no 
use when n, m are > even for a small value of s (sNa 2 — 0.06). The results are equally bad for the other models and 
are not shown here. 

Consider now the case of a realistic model of Na2 18] where both the upper and lower potentials are properly treated; 
the b 1 3 II U to l 3 II g transition has been chosen as an example. The results using the classical Wigner approximation 
(solid line) and the uniform semiclassical approach (dashed line) are compared in Fig. 6(a)-(f). Here the classical 
calculations were carried out using Eq. I|15l) with realistic Na2 potentials, whereas the semiclassical Franck-Condon 
factors were obtained by numerically evaluating 

(E-\j)= f dq(E-\q)(q\j), j=n,m, (34) 

using Gauss-Legendre quadrature. Here, q is the Na-Na separation, (q\j) are the bound vibrational wave functions of 
the initial electronic state calculated quantum mechanically using the Renormalized Numerov method |19|. and (q\E~) 
are the continuum wave functions of the excited potential energy surface, calculated using the uniform semiclassical 
approximation |20j,[2l| 



A O) 2 



(35) 
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FIG. 4: Phase space densities ^5,4 (dashed line) and 5(E — H)w (solid line) at fixed p 
systems. Vertical line corresponds to the classical limit S(q — q(p, E)). 



and E = Eo for three different 
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FIG. 5: Comparison between the classical Wigner and exact quantum (solid line) and quantum corrected Wigner (dot-dashed 
line) results for the n — 5,m = 4, /i = /iNa 2 and ui — 158.91 cm -1 . 



where 



dq'K(q') 



2/3 



(36) 



with 



K(q)=p(q)/h, p(q) = {2 f x(E-V(q))} 1 *, 



(37) 
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5 = 



poo 

l [K(q) - k]dq- ka + -, k = lim K(q), 

Ja 4 



(38) 



and a is the classical turning point, which satisfies V(q = a) = E. As can be seen in Fig 6 the results clearly indicate 
that the classical Wigner representation works very well for the low-lying levels [n,m < 5) of the initial electronic 
state. The results for very high n and m are in poorer agreement and are not presented here. This is a direct 
consequence of the use of the harmonic approximation for the ground state. Specifically, for n, m > 5, anharmonic 
corrections to the ground state should be included in order to account for the derealization of the ground state 
wavefunctions 221. 
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FIG. 6: Comparison between the classical Wigner (solid line) and the uniform semiclassical approximation (dashed line) 
calculations carried out for the realistic Na2 potentials: (a) n = m = 0; (b) n — m = 2; (c) n — m = 4; (d) n = l,m = 0; (e) 
n — 2, m — 0; (f ) n = 5, m = 4. 



V. MULTI-PRODUCT ARRANGEMENT CHANNELS 



Consider now the multi-arrangement channel problem, i.e., the case where photodissociation results in the formation 
of two different chemical products, e.g., A + BC <— ABC — > AB + C. In this case our main focus is on obtaining 
cross sections into specific channels. 

Quantum mechanically the Hilbert space of a typical multi-arrangement channel scattering problem can be parti- 
tioned as follows 1141: 



B. 



(39) 



Here, / is the unit operator, R r projects onto states that correlate asymptotically with all states in channel r, and B 
projects onto bound states. This allows the total cr n rn (E) to be written as 



a n ,m(E) = J2 TT [RrS(E~H)\n)(r 

r 



Tr 



B5(E - H)\n)(m\ 



(40) 
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The channel-specific cross section of interest in this section is given by Eq. J7J in the Franck-Condon approximation, 



i.e. 



= H%T±[Rr6(E -H)\n){m\ 



(41) 



Classically, these operators R r and B correspond to various types of classical trajectories that occur in photodis- 
sociation: trajectories that start in the region of the excited polyatomic (upon excitation) and dissociate to the r 
channels, and those that do not dissociate and remain bound. 

Equation (|41|l . in the Wigner representation, assumes the form 



T {r) 



(E) 



Tr[{R r 5(E-H)} w [\n)(m\} w 
dz[R r 6(E - H)] w (z)p n ,m{z) 



(42) 



where z denotes all coordinates and momenta, q,p. 

To evaluate the channel-specific cross section requires that we approximate the term [R r 5(E — H)]w- In general, 
the Wigner transform of the product of two operators admits a small U expansion |23| 



(AB) w (z) = A w (z)B w {z) + y {A w ,B w } p (z) + 0{h 2 ), 
where {•, •} is the Poisson bracket. Therefore, the channel specific cross section can be approximated by 



^n]mi E ) ~ / dzp, hm (z) 

x [R r ] w (z)S(E - H) w (z) + y {[R r ] w ,5{E - H)w} p (z) 
Equation 144|) can be rewritten by employing the cyclic invariance of the trace, so that 



(43) 



(44) 



(E) 



dz [R r ]w{z) 

5(E - H) W (z)p ntm (z) + y ^S(E - H)w,Prum} (z) 



(45) 



This form is more natural since the term that selects the channel, [i? r ]vK(z)) acts in the same manner on both terms 
in the integral. 

We can implement this channel selection as follows: we consider the trajectory that emanates from an initial point 
z; if the trajectory ends in channel r it contributes to an)m{E). Alternatively, it contributes to channel r' ^ r, or is 

(r) 

bound, both of which are ignored in the a n ' m computation. 

Both the magnitude and phase of <Jn}m (E) are important in coherent control. Hence we note that this term can have 
an imaginary part if n ^ m, since the integration over momentum is now constrained and arguments based on the 
odd vs. even nature of the integrand do not apply. To see this more clearly, consider <Jn}m{E) = {m\R r 8(E — H)\n). 
The wavefunctions \ j) are real, and therefore if On}m{E) is to have an imaginary component, the operator R r S(E — H) 
must be non-Hermitian. Since each of R r |14| and 8{E — H) are individually Hermitian, the operator R r 8(E — H) is 
non-Hermitian only if [R r , S(E — H)] ^ 0. This is indeed the case, as can be seen by taking the matrix element with 
respect to eigenfunctions of the Hamiltonian H, to find that 



(E \[Rr,S(E- H)]\E ) = {E \R r 5(E — H)\E") - (E\S(E - H)R r \E ' 
= (E'\R r \E")(6(E - E) - 5{E - E)) ^ 0. 



(46) 



However, even though Un}m{E) can have an imaginary component, the sum of the imaginary contributions over 
all the channels must remain zero since u n ,m{E) is real. For a two channel system, this implies that any imaginary 
contribution in channel 1 must be equal in magnitude and opposite in sign to the imaginary contribution from channel 
2. We further note that the diagonal <Jn}n(E) is 



ai r UE) = (n\R r S(E-H)\n), 
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k 

= 2|(n|£?,k,r-)| a 
k 

and is therefore real. 

Consider now the second term in Eq. . The derivatives of the delta function introduced via the Poisson bracket 
result in the rapid oscillation of the integrand which, as we have verified numerically, yields a final contribution to the 
integral that is essentially zero. Therefore, we set the Poisson bracket term to zero, and the channel specific (Jn}m{E) 
becomes 

<r$n{E) - Jdz [R r ] w (z)S{E - H) w {z)p n , m {z). (47) 

Equation (|47|l has the form of two overlapping phase space densities. Specifically, the density p n ,m(z) overlaps 
[R r ]wS{E — H)w, where the latter is the phase space density of states on the energy shell at energy E that decays 
to product in channel r. To successfully approximate Eq. l|T7|) requires a classical approximation to the Wigner 
transform [R r ]w[&{E — H)]w, discussed below. 



A. Computational Results 

The lowest dimensionality problem of this kind, useful for examining the utility of the semiclassical approximation 
under consideration, is the collinear photodissociation of ABC as A + BC <— ABC — ► AB + C, where each of A, 
B, C denote an atom or molecular fragment. We consider this problem below, where the electronic ground state is 
assumed to be well approximated by an harmonic oscillator. 

For the two degrees of freedom case we use the notation z = (pi,p2, qi, 32), P = (pi^Pz), Q = ((71,(72), n= {nx,ri2), 
m = (mi, 7712), where Pi, qi denote the momenta and coordinates of the two degrees of freedom. The two dimensional 
harmonic oscillator initial vibrational state is given by 

(qiK)((72M - N ni H ni {^r 1 q 1 )e- a ^' 2 N n2 H n2 {^r 2 q 2 )e- a ^' 2 1 (48) 
and the two dimensional Wigner function can be written as the product of two one dimensional Wigner functions 

PnM z ) = Pni,n a ,mi,m 2 (2) = P^umi (Pl> Ql)Pn^ fe, 32). (49) 

The excited state Hamiltonian is 

In the computations below it proves advantageous to numerically approximate the delta function as 

1 ( -(E - H(z)) 2 ) , s 

S(E-H) W ~5(E-H(z))~^=cxpl-^ 4e J ' (51) 

where e is chosen small. The final form for the approximate channel specific term in two dimensions then becomes 

1 f ( — (E — H (z)) 2 } 
{E) ~ J dz [Rr(z)}wPn u n 2 , mi , m2 (z)exp<^ — 1 . (52) 

Results were also obtained by expanding the delta function to reduce the dimensionality of the integrand from four 
to three: 

^-^ )) = 1^ + T^ (53) 

\pi/M 2 \ \PilM 2 \ 

where 



Pt = pf(qi,q 2 ,Pi,E) - ±J2M 2 [E - p 2 /(2M 1 ) - V(q u q 2 )]. (54) 
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so that the channel specific term is now given by 

a ni n 2 mi m 2 ( E ) = [ d 1l dq 2 dpi [R r (qi , q 2 , Pi , P2 )] WPm .n 2 ,mi ,m 2 (<7l , <?2 , Pi , P2 ) , (55) 

p 2 \P2\ 

where the integral over T requires that E > V(qi 1 q 2 ) + p\/{2Mi), and the sum is over p 2 = pf- 

Monte Carlo integration of Eq. i|52[l and (|55|l showed that Eq. I|52|) converged with fewer trajectories, even though 
the integral is of higher dimension. Further, the ability to smooth the integral by increasing e allows qualitative 
estimates of the form of the channel specific cross sections with a small number of trajectories. Some sample results 
are provided below. 



VI. APPLICATION TO C0 2 AND CH 2 BRI 

The method was first applied to the photodissociation of collinear CC^- The coordinates are denoted R — ro-c, r = 
rc-o- The system is initially in the bound state x E'&(r, R), with equilibrium separation x s — R = f = 2.20 a.u. The 
ground potential surface is harmonic in the normal mode coordinates x s and x a |24|. with parameters: 

x s = (R + r)/2; x a = (R- r)/(2 7 ); 7 = [1 + m c /(2m )], (56) 
M s = 2m ; M a = m c (l + m c /(2m )), (57) 
a s = u s M s /H; a a = uj a M a /h. (58) 

The two product channels are denoted OC+O (channel 1) and O+CO (channel 2). The coordinates qi = x s , q 2 = x a 
are related to the bond length coordinates by 

r-f = qi - 752, 
R-R = qi+iq 2 . (59) 

The multidimensional integrals in Eqs. (|52|l and (|55|l for both this case as well as the case of CH^IBr discussed below 
were carried out for all systems below using Monte Carlo box Muller transformation [23] and 10 6 trajectories. The 
value of e was chosen as e = 5 x 10~ 6 for the case of CO2 and e = lx 10~ 8 for CH^IBr. 

Results for CO2 are shown in Figs. 7 and 8. Figure 7 compares the total cross section ctoooo(^) calculated using 
Eq. (|52H . computed with R r (z) = 1, to results of a time dependent formalism utilizing the stationary phase Herman 
Kluk (SPHK) propagator j25j. Trajectories in the SPHK procedure were followed only long enough to capture the 
initial dispersion from the Franck-Condon region, giving the direct part of the cross section. Although the classical- 
Wigner result is seen to be shifted slightly from the time dependent result, the method is seen to display the essential 
features of the direct part of the photodissociation cross section. An analagous calculation by Eckhardt and HiipperQ 
for water produced a result similar to that shown in Fig. 7. 




1 ■ ■ ■ i 

1.5 2 2.5 3 3.5 

Energy (a.u.) 



FIG. 7: CO2 results for aoooo(E) (total result) using the time dependent (autocorrelation) formalism (dash) and the time 
independent (solid) formalism. 

Figure 8 shows a^Q 00 (E) , ajf^^E) , and cr^ 00 (E) for CO2. In this case the diagonal cross sections are strictly real. 

The off-diagonal cross section Cqioo^) d° contain an imaginary component in both channels, but they are seen to be 
equal and opposite in magnitude, verified by computing the total cross section which is again strictly real. 
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FIG. 8: CO2 results. The solid line is the sum from both channels, the dash lines are the two channel contributions, (a) 
Re[<r« (£)] = 4» (£); (b) Re[a^ 01 (E)] = a^ 01 (E); (c) Re[a^ 00 (E)]; (d) lm[a^ 00 (E)]. 



To test the utility of this approach on a system where the two channels are dissimilar, we consider CH^Brl where 
channel 1 is CH2B1- 



I and channel 2 is CH2I 

T^(rCH 2 -Br, ?*CH 2 -l) = 



- Br. The upper potential energy surface is given by [26 
[VKt-ch 2 -i) + V^BrCrcHi-B,)]/^) + 

[^Br(r CH2 -Br) + (r C H 2 -l)] (1 - f(x)) 



(60) 



where 



T/Morsc 
^CH 



Brl r CH 2 -Br 



^CT 2 -Br ex P{-aBr(V C H 2 -Br - ?CH 2 -Br)} 

x[exp{-a B r(r C H 2 -Br - ?CH 2 -Br)} - 2 ] 



^c M h° 2 -i(^h 2 -i) = ^^expj-ajCrcH.-i 



' CH 2 



-i)} 



Vi(rcH 2 -i) 

VBr(7CH 2 -Br) 



x [exp {-ai(rcH 2 -i - r e cn2 
Aiexp{-/3ir C H 2 -i} , 

A Bl exp {-/3 Bl TCH 2 -Br} ) 

1 



-i)}-2] 



(61) 

(62) 
(63) 
(64) 

(65) 
(66) 

3.6850 a.u., 

A Bl = 0.27, (3 Bl = 0.35, £>ch 2 -i = 0.0874 a.u., aj = 0.87094 (a.u.)" 1 , r^.j = 4.04326 a.u., Aj = 0.37, fa = 0.3. 

The initial state is taken to be harmonic in normal mode coordinates. These coordinates, qi,q2, are related to the 
bond length coordinates rcH 2 -Br an d ^ch 2 -i by 



1 



exp {a(x 

rCH 2 -Br 



0-5)} : 



f"CH 2 -I + ?X:H 2 -Br 



The parameters for this surface are: a = 30, ^ch -Br = 0-1069 a.u., as r = 0.9154 (a.u. 



CH 2 -Br 



^CH 2 -Br - »"CH 2 -Br = CnQi + C 2 iq 2 
?*CH 2 -I - ^CH 2 -I = Ci2<7l + C22 92 

where c u = 0.552747; c 2 i = 1.09614; c 12 = 0.788417; c 22 = -1.03 8 93; r CH2 - B r = 
parameters for this system are: 

Mi = 162614.16 a.u.; M 2 = 27238.15 a.u., 

c*i = 149.3978 a.u.; a 2 = 98.8412 a.u., 

flux = 201.638 cm" 1 ; Tiu> 2 = 796.43112 cm" 1 . 



r C H 2 -i = % 2 _i- 



(67) 
(68) 

The 



(69) 
(70) 
(71) 
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Results for several cross sections and interference terms are shown in Figs. 9 to 11. The diagonal cross sections 
are strictly real. For this system, unlike CO2, the real part of the channel specific results can be different in the two 
channels. 
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FIG. 9: CHbBrI results. The solid line is the sum from both channels, the dash lines are the two channel contributions, (a) 
M4 k oUE)} = ^ooo(-B); (b) Re[*%(E)] = a{ k \ (E); (c) Re[a^ (E)]; (d) Ma^ 00 (E)]. 




FIG. 10: CHhBrI results, comparing the Wigner-classical method (solid) to the quantum (dash) for the CtFBr + I channel, 
(a) Re[<T$ 00 (E)) = <t$ Q0 (E); (b) Re[a^ 10 (E)} = a$ 00 (E); (c) Re[<r§ 00 (E)}; (d) lm[a^ 00 (E)}. 

Figure 9 shows the structure of the channel contributions to various cross sections and interference terms. Note 
that, as required, for the case of the off-diagonal cross section cr^oiE) is interesting (cf. Fig. 9 ), even though the 
two channels produce different real parts, the imaginary parts in the two channels again sum to zero. The deviation 
from zero is very small, indicating reliable convergence of the Monte Carlo sums. 

Figures 10 and 11 show a comparison of the classical- Wigner method with a full quantum mechanical calculation. 
The structure of the cross section in the channels is seen to be well reproduced the classical- Wigner method, providing 
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FIG. 11: CHhBrI results, comparing the Wigner-classical method (solid) to the quantum (dash) for the CH2I + Br channel, 
(a) Re^ (£)] = Cqooo (E); (b) Re[a^ 10 (E)] = a{ 2 [ (E); (c) Re^f^^iJ)]; (d) lm[a[l\ )0 (E)]. 

support for the conclusion that this approach gives reliable results for both the real and imaginary contribution. 
Interestingly, the latter arises via the non-Hcrmitian character of R r S(E — H). 

Finally, Fig. 12 shows an example of the utility of larger values of e in Eq. Q51f) . Specifically, using larger values of 
e allows for a qualitative estimate of the desired integrals using far fewer trajectories. 
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FIG. 12: CH 2 BrI results for Cqooo (E) (CH 2 Br + I) with different e. Dash: e = 5 x 10~ 7 , 2 x 10 4 trajectories; Solid: e = 10" 
10 6 trajectories per energy. 



VII. SUMMARY 



We have considered a Wigner-based classical approximation [Eq. I|14|) ] to compute the terms 
^2 r {E, k~\[ifi\n)(m\fj,fi\Ei fc~), with n 7^ m, that are central to the interference contributions characteristic of co- 
herent control. In the case of a single open product arrangement channel the accuracy of this formula for the 
Franck-Condon transitions onto a linear potential, along with the dependence of the method on parameters such as 
system reduced mass, slope of the upper potential energy surface, and vibrational frequency of the lower electronic 
state were examined. The results were found to be in excellent agreement when the newly introduced parameter 
s n ,m [Eq. (|33() ] was less than unity. A comparison of the classical and uniform semiclassical approximations for the 
transitions between the realistic potential energy surfaces of Na2 demonstrates that higher values of m, n require use 
of the anharmonic Wigner function for the ground state, an effect not explored in this paper. 
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This approach was also applied to the mult i- channel problem where interference terms are complex. In this case, the 
result is still the overlap of two time-independent phase space densities, but determining the density associated with 
a particular channel required evaluation using classical trajectories. The non-Hermitian character introduced by the 
projection operator R r onto channel r allowed for the successful reproduction of the entire complex term. We regard 
it as particularly encouraging that use of classical trajectories in conjunction with the complex Wigner transform of 
the term suffices to produce the imaginary part of the interference contribution with such accuracy. We note, 

however, that the method is not expected to produce resonance structures in the cross section, being most useful for 
systems where the dynamics is direct and short-lived|8l l9l lllf. 

Applications to full 3 dimensional photodissociation computations and to scattering (for bimolecular coherent 
control are underway. 
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